Long-term prognostic impact of paravalvular leakage on coronary artery disease requires patient-specific quantification of hemodynamics

Transcatheter aortic valve replacement (TAVR) is a frequently used minimally invasive intervention for patient with aortic stenosis across a broad risk spectrum. While coronary artery disease (CAD) is present in approximately half of TAVR candidates, correlation of post-TAVR complications such as paravalvular leakage (PVL) or misalignment with CAD are not fully understood. For this purpose, we developed a multiscale computational framework based on a patient-specific lumped-parameter algorithm and a 3-D strongly-coupled fluid–structure interaction model to quantify metrics of global circulatory function, metrics of global cardiac function and local cardiac fluid dynamics in 6 patients. Based on our findings, PVL limits the benefits of TAVR and restricts coronary perfusion due to the lack of sufficient coronary blood flow during diastole phase (e.g., maximum coronary flow rate reduced by 21.73%, 21.43% and 21.43% in the left anterior descending (LAD), left circumflex (LCX) and right coronary artery (RCA) respectively (N = 6)). Moreover, PVL may increase the LV load (e.g., LV load increased by 17.57% (N = 6)) and decrease the coronary wall shear stress (e.g., maximum wall shear stress reduced by 20.62%, 21.92%, 22.28% and 25.66% in the left main coronary artery (LMCA), left anterior descending (LAD), left circumflex (LCX) and right coronary artery (RCA) respectively (N = 6)), which could promote atherosclerosis development through loss of the physiological flow-oriented alignment of endothelial cells. This study demonstrated that a rigorously developed personalized image-based computational framework can provide vital insights into underlying mechanics of TAVR and CAD interactions and assist in treatment planning and patient risk stratification in patients.

Patient-specific lumped parameter algorithm for complex valvular, vascular, mini-vascular and ventricular diseases. We have previously developed a non-invasive diagnostic computational-mechanics framework for complex valvular, vascular and ventricular disease (called C3V-LPM for simplicity), described in detail elsewhere 21 . In this study, we further developed the C3V-LPM to enable the quantification of local and global hemodynamics in patients with mixed and complex valvular, vascular, mini-vascular and ventricular diseases (known as C3VM-LPM) (Fig. 1, Table 2). The developed framework uses limited input parameters that can all be reliably measured non-invasively using DE, CT and a sphygmomanometer. The previously created model, C3V-LPM, was validated against clinical catheterization data in forty-nine AS patients with a substantial interand intra-patient variability with a wide range of disease 21 . In addition, some of the sub-models of the patientspecific lumped parameter algorithm have been used previously 12,14,25-31, , with validation against in vivo cardiac catheterization (N = 34) 32,33 in patients with vascular diseases, in vivo MRI data (N = 57) 34 in patients with AS, and in vivo MRI data (N = 23) 35,36 in patients with coarctation and mixed valvular diseases.
Left ventricle. Time-varying elastance, E(t), is a common method to simulate the ventricle muscle stiffness which relates the LV pressure and volume: where P LV (t) , V (t) and V 0 represent the LV time-varying pressure, time-varying volume and unloaded volume, respectively. This elastance function is often represented by the double Hill function, initially proposed by Stergiopulos et al. 38 . This mathematical function is capable of capturing the contraction and relaxation dynamics of the ventricle where τ 1 , τ 2 , m 1 , m 2 , E max and E min are ascending time translation, descending time translation, ascending gradient, descending gradient, maximum elastance and minimum elastance, respectively (see Table 2). The first term (1) Table 2. www.nature.com/scientificreports/ in brackets represents the contraction of the chamber and the second term in brackets represents the relaxation of the chamber. As outlined in Table 2: (1) τ 1 and τ 2 are functions of the cardiac cycle duration (T) and vary for each patient; (2) m 1 , m 2 are constant for all patients (Stergiopulos et al. 38 , Mynard et al. 39 , Seemann et al. 45 ). Parameter values used in the elastance function were adapted to generate physiological waveforms for pressure, volume and flow that can be found in Table 2 38,[46][47][48][49][50][51][52][53][54] .
Left atrium. The same time-varying elastance function, E(t), used for the LV was applied to the LA for the coupling of LA pressure and volume. The parameter values used for the LA are found in Table 2, while the elastance function is defined by Eqs. (2) and (3) 21 .  Table 2. Simulation domain and FSI modeling. Patient-specific LPM simulating the function of the left side of the heart and coronary arteries provided the patient-specific boundary conditions of the inlet and outlets. Mitral valve. The mitral valve (MV) was modeled using a similar technique to the aortic valve which accounts for the net pressure gradient ( PG net | MV ) across the MV during LA ejection. PG net | MV is a function of ρ , Q MV , EOA MV and M MV , which represent the density of the fluid, the transvalvular flow rate, effective orifice area and inertance, respectively.  where EOA| MR is the MR effective orifice area. Figure 7. Grid convergence analysis. (a) TAWSS at the PVL location: TAWSS comparisons for different mesh resolutions at the peak of filling phase (maximum TAWSS): difference in peak TAWSS between mesh#3 and mesh#4 is 3% and between mesh#4 and #5 is 1.03%; (b) Von Mises stress: Von mises stress comparisons for different mesh resolutions of left coronary artery at the peak of filling phase: difference in maximum stress between mesh#3 and mesh#4 is 5.25% and between mesh#4 and #5 is 1.43%; (c) Plots for Maximum TAWSS and Maximum Von Mises stress versus the number of mesh elements. where Q MPV , t ee and T are the mean flow rate of the pulmonary valve, end-ejection time, and cardiac cycle time, respectively. It is important to note the only input flow condition that can be reliably measured using DE in this study is the forward left ventricular outflow tract stroke volume (Forward LVOT-SV). The lumped-parameter model could reproduce this DE-measured metric when Q MPV , the mean flow rate of the pulmonary valve, was optimized.
Coronary arteries. C3VM-LPM was developed to simulate blood flow rate at the outlets of the left anterior descending (LAD) artery, left circumflex (LCX) artery and right coronary artery (RCA), in addition to other regions of the heart and the cardiovascular system. The 5-element electrical circuit used in this study to model each coronary branch was developed by Kim et al. 40 and has been used extensively to generate boundary conditions for higher order coronary models 41,42,44,[55][56][57][58] . The circuits comprised of 3 resistors ( R cor,p , R cor,m , R cor,d ) , 2 capacitors ( C cor,p , C cor,m ) and an embedded pressure source ( P im ) (Figs. 1 and 2). While inductors are often included in the heart and circulatory models, they were not included in this coronary model since the inertial phenomena is not significant in the coronary artery branches 59 . By including an embedded pressure source, this layout has been shown to capture the biphasic nature of coronary flow, in which peak blood flow occurs during the diastole phase rather than during the systole phase 40,59 .
The following ordinary differential equations (ODEs) are generated to model the coronary lumped parameter model 55 :  www.nature.com/scientificreports/ where q in , P in , q out and P out are the blood flow and pressure into and out of the coronary branch. R cor,p , R cor,m , R cor,d are the proximal, medial, and distal resistors while C cor,p , C cor,m are the proximal and medial capacitors. P p , P m and P im are the proximal, medial and intramyocardial pressures. P im is set to be either the left ventricle (LV) or right ventricle (RV) pressure, depending on the coronary artery that it is coupled to. In this study, we used the LV pressure for the left branches (LAD and LCX) and 0.5P LV 40 to create the RV pressure for the right branch (RCA).
Determining arterial resistance and compliance in coronaries. Total coronary resistance The total coronary resistance was derived from a relationship between blood pressure and blood flow, where the mean flow rate to the coronary arteries was assumed to be 4.0% of the cardiac output 40,41 : (11) q in = P in − P p R cor,p (12) q in = C cor,p dP p dt + q m (13) P p = q m R cor,m + P m (14) q m = q out + C cor,m dP im dt (15) P m = q out R cor,d + P out    www.nature.com/scientificreports/ Coronary vessel resistance and compliance The total coronary resistance was divided between each of the branches based on a generalization of Murray's law 43 , which relates resistance to vessel diameter: R cor,total where j = {LAD, LCX or RCA} www.nature.com/scientificreports/ where R cor,j is the total coronary resistance in the desired vessel and A i is the cross sectional area of each of the coronary vessels 40 . Further division of the total vessel resistance into the 3 resistive elements in the circuit was based on the work of Sankaran et al. 42 :   www.nature.com/scientificreports/ where C cor,j is the left coronary vessel compliance, C L cor,total is the total left coronary compliance and A i is the cross sectional area of each of the left coronary branches 40 . A manual tuning process was utilized to determine total left coronary compliance value that lead to physiological coronary flow waveforms 44,61,62 .
The compliances were then divided across the 2 capacitors based on the following relationship 42 : where C cor,j,p and C cor,j,m are the proximal and medial capacitors. The same process was employed for the right coronary vessels.
Input parameters. The following patient-specific parameters were inputs for the C3VM-LPM algorithm: forward left ventricular outflow tract stroke volume (Forward LVOT-SV), cardiac cycle time, ejection time, aortic valve effective orifice area (EOA), mitral valve EOA, ascending aorta cross sectional area, left ventricle outflow tract area, EOA during aortic regurgitation and EOA during mitral regurgitation measured by DE as well as brachial systolic and diastolic pressures measured by a sphygmomanometer. All the details about patient-specific parameter estimation were described in 21 . In addition, coronary geometry dimensions (left main coronary artery (LMCA) average diameter, LAD coronary artery average diameter, LCX coronary artery average diameter (21) C cor,j,p = (0.11)C cor,j C cor,j,m = (0.89)C cor,j www.nature.com/scientificreports/ and right coronary artery (RCA) average diameter), measured from the patient-specific reconstructed coronary artery geometry using CT data, are the input parameters for the C3VM-LPM algorithm.
Computational algorithm. MathWorks Simscape (MathWorks, Inc.) was used to formulate and solve the system of ordinary differential equations (ODEs) which govern the lumped parameter circuit. Additional functions were written in Matlab and Simscape to supplement and enhance the Simscape code. The ode23t trapezoid rule variable step solver with an initial time step of 0.1 ms was used to solve these ODEs. Initially the voltages and currents in the capacitors and inductors were set to zero and the model was run for ~ 150 cycles to reach a steady state. For the patient specific optimization, the residual criterion was set to 10 −6 .
We generated a signal to model LV elastance using a double Hill function representation of a normalized elastance curve for human adults 21 . The LV pressure, P LV , calculated using the initial values of the model input parameters from Table 2, and the time-varying elastance (Eq. 1), were used to compute the instantaneous LV volume, V(t). Subsequently, the time derivative of the instantaneous LV volume was calculated to find the LV flow rate. This approach was also applied to obtain the volume of the left-atrium, pressure, and flow rate.
Patient-specific response optimization. To correctly patient-variability, four parameters of the model were optimized such that the lumped-parameter model reproduced the physiological measurements obtained from the patient. Simulink Design Optimization toolbox was used to optimize the response of the lumped-parameter model based on Matlab's fmincon function.
Since Forward LVOT-SV can be measured reliably using DE, Q MPV was tuned to minimize the difference between the Forward LVOT-SV calculated by the model and the one measured using DE in each patient: where D LVOT , A LVOT , and VTI LVOT are LVOT diameter, LVOT area, and LVOT velocity-time integral, respectively 21 .
In the second step, R SA , C SAC , and C ao were optimized so that maximum and minimum aortic pressures were respectively equal to the systolic and diastolic pressures measured using a sphygmomanometer in each patient. Since the left ventricle experiences the total systemic resistance rather than the individual resistances, and the www.nature.com/scientificreports/ systemic arteries resistance, R SA , is one order of magnitude greater than both the aortic resistance, R ao , and systemic vein resistance, R SV , for the sake of simplicity we considered R ao and R SV as constants and optimized R SA as the main contributor of the total systemic resistance. C ao was considered to be 0.6 of C SAC because 60% of the total arterial compliance resides in the proximal aorta 63 . End systolic volume (ESV) or end diastolic volume (EDV) measured by DE was fed to the lumped-parameter model to adjust only starting and ending volumes in the P-V loop diagram. For this purpose, the Biplane Ellipsoid model was used to calculate the instantaneous LV volume at the end of diastole or systole as follows 21 : where A 1 , A 2 , L 1 and L 2 are LV area measured in the apical four-chamber view, LV area measured in the apical two-chamber view, LV length measured in the apical four-chamber view and LV length measured in the apical two-chamber view, respectively.
In addition, we conducted an extensive parameter sensitivity analysis that revealed changes in pulmonary parameters (e.g., C PVC ) have negligible effects on the model output variables. Therefore, we did not include these pulmonary parameters in the parameter-optimization process and considered them as constants given in Table 2.
Fluid-structure interaction simulation study. The blood flow inside the coronary arteries was simulated using similar 3-D FSI set-up as described in our previous works 18,19 , using open-source FOAM-Extend Table 1. Baseline and post-TAVR patient characteristics. C3VD patients (n = 6, mean ± SD) Pre-TAVR (n = 6, mean ± SD) Post-TAVR (n = 6, mean ± SD) www.nature.com/scientificreports/ library 64 . The transcatheter aortic valve (TAV) frame and aortic wall are assumed to be rigid during diastole 65,66 . All details about governing equations, FSI method and modeling were presented in Supplementary Materials.
Boundary conditions. Previous studies 65,67,68 have used time-dependent pressure waveform of ascending aorta and ventricle (assuming fixed and rigid valve and aorta) to obtain the main hemodynamic features of PVL during diastole. However, the accurate and patient specific time-dependent pressure waveform is necessary for any CFD or FSI simulation as the whole topology of fluid domain is affected by pressure waveforms, and that's why a patient-specific pressure waveform is crucial for an accurate simulation. Our patients specific lumped parameter algorithm generated boundary conditions (Figs. 1 and 2) non-invasively to provide 21 : (1) the time-dependent pressure waveform of ascending aorta during diastole which was applied as inlet boundary condition; (2) the time-dependent coronary flow waveforms which were applied as coronary outlet boundary condition; (3) for PVL, the ventricle pressure was applied as outlet at the leakage area location (leakage area was measured and located based on short axis DE, Fig. 3b) to provide patient-specific PVL pressure gradient (pressure difference between ascending aorta and ventricle during diastole). For the PVL simulation, we used the same approach in previous studies 67,68 . The leakage surface is the interface wall between aortic root (behind the TAV stent) and LVOT which is considered for outlet boundary condition. This surface is in the vicinity of the outer region of TAV stent. Therefore, the blood flow was not free (i.e., pressure zero in the outlet) to move towards ventricle. Instead, the blood flow was driven towards the LVOT (and ventricle) by pressure gradient between ascending aorta and left ventricle during diastole (Fig. 4). All the simulations were performed during diastole and the aortic valve was therefore assumed to be rigidly closed since the large deformation of valve and aorta occurs mainly during systole, while they remain relatively motionless during diastole with negligible effect of detailed closing shape of the leaflets on the PVL 68,69 . It is important to note that PVL occurs only during diastole 70 and previous studies 67,68,71 have performed the simulation only for diastole under assumption of rigid valves/aorta 65,67,71 and found that this assumption does not affect the conclusions of their studies. They have also validated their results with experimental and clinical data 65,67,71 . We followed the same approach 65,67,71 in our study and performed validation ("Validation: Doppler-based LPM and FSI framework versus clinical Doppler echocardiography data" section) with clinical DE velocity magnitude. However, we acknowledge that this a limitation of our study and we addressed this in the limitation section. Yet, it is worth mentioning that in the absence of essential characterization of patient-specific material properties required for FSI simulation of valve flow (which is the case in all ongoing FSI simulations 69,72,73 ), the results of FSI simulation could be an incorrect representative of the flow 66 . Therefore, whether the FSI simulation of the valve and aorta with such limitations improves the results is still debatable especially if the end goal is to provide a patient-specific framework 74 . Figure 8a to f compares the peak PVL velocity simulated using our computational framework and DE data for two patients as a sample (8a and 8d: regurgitant flow waveform 8b: parasternal short axis view of PVL jet; 8e: parasternal long axis view of PVL jet). The simulated peak velocities are in a good agreement with the ones measured by DE in both patients with a maximum relative error of 8% for the peak velocity at the beginning of diastole phase (early filling). For the whole diastole phase, the results show good agreements between velocity calculated using the computational framework and the ones measured using DE in both cases investigated in this study.

Effect of anatomic and deployment characteristics on aortic root and neo-sinus local hemodynamics (post-TAVR).
The blood flow vortical structure and stagnation in the aortic root, sinus of Valsalva and neo-sinus region depends on the ascending aorta and ventricular pressures, aortic root geometry, aortic valve geometry, stent height, deployment angle and coronary ostium location. We investigated hemodynamic metrics computed by our computational framework (Figs. 9 and 10) as follows: Vortical structure. It has been shown that for a TAVR without PVL, the coronary flow influences the flow patterns of aortic root and neo-sinus and favors the transfer of blood flow towards ostium during diastole 75 . However, our results showed that in the presence of PVL, the aortic root vortices will not favor the transfer of blood flow towards ostium in the aortic root and neo-sinus region. As shown in Figs. 9a,b,c and 10a,b,c, for patients #1 and #2, the mainstream of PVL flow originates from ascending aorta towards the leakage orifice behind the stent and between LCC and RCC leaflets with a maximum of 2.05 m/s and 3.22 m/s for patients #1 and #2 respectively. However, the maximum velocity between left ostium and stent was 1.53 m/s and 0.42 m/s for patients #1 and #2 respectively. This can be explained by the fact that the size of the gap between the edge of stent frame and the ostium is smaller for patient #1 than patient #2, leading to higher divergent velocity towards the leakage area. For both patients, a vortex forms in the neo-sinus region of all the leaflets (LCC, RCC and NCC) as shown in Figs. 9a,b,c and 10a,b,c. Our results showed this vortex is very different for LCC, RCC and NCC and for different patients. For patient #1, the vortex arises from the leaflet surface at early diastole and dominates the whole neo-sinus region at late diastole, leading to an efficient washout of blood flow from the LCC (Fig. 9a). However, during the whole diastole, the center of vortex remains close to the stent edge for RCC (Fig. 9b), and for NCC (Fig. 9c), the center of vortex remains close to the upper commissure and vanishes at late diastole. For both RCC and NCC (Fig. 9b,c), the vortex does not move down to reach the leaflet surface, leading to a reduced washout of blood flow. For patient #2 though, the vortex center remains distant from the leaflets for LCC (Fig. 10a) and NCC (Fig. 10c) at early diastole and gets closer to the leaflets for RCC (Fig. 10b). Although in mid diastole, the www.nature.com/scientificreports/ vortex size in LCC increases, and does not dominate the whole neo-sins similar to the vortex for patient #1. In other words, for patient #2, vortices aid the washout in RCC more than LCC and NCC.

Stagnant and low-velocity flow.
For patients #1 and #2 as a sample, NCC neo-sinus had higher regions of stagnant flow than RCC and LCC neo-sinuses; 0.24 cm 3 at early diastole and 0.45 cm 3 at late diastole for patient#1 (Fig. 9d), and 0.55 cm 3 at early diastole and 0.903 cm 3 at late diastole for patient#2 (Fig. 10d). For patient #1 the RCC had significantly higher stagnant flow than LCC for the whole diastole; 0.23 cm 3 and 0.042 cm 3 at early diastole and 0.12 cm 3 and 0.064 cm 3 at late diastole for RCC and LCC respectively. However, for patient #2, the LCC had slightly higher stagnant flow at early diastole and RCC had slightly higher stagnant flow at late diastole; 0.075 cm 3 and 0.067 cm 3 at early diastole and 0.141 cm 3 and 0.153 cm 3 at late diastole for LCC and RCC respectively. Interestingly, although patient #1 had more severe PVL than patient #2, the stagnant flow volume in the neo-sinus region was almost 2 folds larger for patient #2 than patient #1 (LCC and NCC). This can be explained by the fact that the blood stasis depends not only on the PVL severity, but also on the patient-specific aortic root geometry, ascending aorta and LV pressures and the deployment details of TAVR. In other words, our results showed that PVL severity alone cannot reveal the risk of thrombosis in the neo-sinus region.
Aortic root wall shear stress. Wall shear stress, as a tangential force induced by blood flow, has a major influence on regulating endothelial function 76 . In general, very high wall shear stress (typically higher than 3 Pa) could contribute to tissue rupture 77 . PVL could disturb the flow in the aortic root sinus after TAVR, leading to increased wall shear stress. As an example, the maximum local TAWSS at the aortic root was increased drastically after TAVR for patients #1 and #2 (Fig. 11); from 0.11 Pa and 0.08 Pa pre-TAVR to 12.6 Pa and 11.8 Pa post-TAVR for patient #1 and patient #2 respectively. Such considerably high TAWSS might be a concern for patients who received TAVR. Moreover, our finding showed that the distribution of wall shear stress at the aortic root is very different for each patient, depending on the characteristics of TAVR deployment and aortic root geometry (Fig. 11).

Coronary arteries blood flow and tissue assessment (pre-TAVR and post-TAVR).
In the presence of PVL after TAVR, the supplied blood flow through the coronary arteries is altered. We investigated the metrics of tissue (solid domain) and flow (fluid domain) computed by strongly coupled FSI algorithm as follows: Coronary arteries von-mises stress. Although there is no cut-off threshold available in the literature for the rupture of arterial wall von-Mises stress, an average stress of 0.3 MPa has been reported to initiate the first crack in the artery wall 78 . As shown in Fig. 12, the distribution of von-Mises stress and its maximum, is different for inner and outer layers of tissue. The maximum von-Mises stress magnitude for all coronary branches in our study was less than 0.3 MPa during diastole for both pre-TAVR and post-TAVR. A universal reduction in maximum von-Mises stress was observed after TAVR for all coronary arteries; As an example, 26.3% reduction for left coronary branches and 11.11% reduction for right coronary in patient #1, and 10% reduction for left coronary branches and 14.3% reduction for right coronary in patient #2.
Coronary arteries wall shear stress. Endothelial cells which are exposed to low wall shear stress display a proinflammatory state, which is associated with plaque progression 76,79 . Although providing an exact cut-off value for low wall shear stress is still challenging, some studies suggested that wall shear stress lower than 1 Pa 76 or 1.2 Pa 79 is associated with higher plaque progression in a further serial study of coronary atherosclerosis. We calculated the wall shear stress over diastole for all patients in pre and post TAVR states. For patients#1 and #2 as examples, local and maximum wall shear stress were decreased for all coronary branches (LCX, LAD and RCA) as shown in Fig. 13. For patient #2, the maximum wall shear stress slightly reduced for LAD and LCX branches; 8.5% and 12.5% for early and late diastole. However, for RCA, maximum wall shear stress decreased from 1.05 to 0.87 Pa (17%) at the peak diastole and from 0.5 to 0. Computed global hemodynamics. Cardiac function. LV workload represents the total energy required by the ventricle to eject blood, and is an effective metric of LV load and clinical state 12,14,18,19 . For patients #1 and #2 for example, despite the reduction of transvalvular pressure gradient, the LV workload increased after TAVR due to the presence of PVL (Fig. 14a); 35% and 18.67% increase in workload after TAVR for patient #1 and #2 respectively. Although the LV pressure decreased post-TAVR, severe PVL contributed to a shift from ventricular pressure overload to a ventricular volume overload.
Circulatory function. Systemic arterial compliance (SAC) is an index for predicting vascular disease states. For patients with AS, a low SAC (lower than 0.64 ml/m 2 /mmHg) is associated with increased risk of morbidity 80 . As shown in Fig. 14a, SAC improved for patients #1 and #2 after TAVR, with SAC increasing to > 1 (ml/mmHg) for both patients after intervention. Increased aortic pressure is expected after TAVR as a result of the removal of AS obstruction 81,82 . As shown in Fig. 14a, maximum aortic  www.nature.com/scientificreports/ aortic pressure increased only 5.1%. Moreover, maximum left atrium pressure reduced by 39% for patient #1, while the change was almost negligible (less than 3% increase) for patient #2.
Coronary circulatory function. Inadequate coronary flowrate and coronary hypoperfusion could lead to exacerbated heart failure 83 . It has been shown that the TAVR deployment characteristics (such as implant depth, angle and PVL) could affect the coronary flow [83][84][85][86] . As shown in Fig. 14b, for all patients in our study, although the perfusion pressure has increased after TAVR, the PVL and flow disturbance in the aortic root significantly reduced the flowrate in almost all coronary branches. For example, maximum flowrate was reduced by 34% and 37% in LAD and RCA branches of patient #1 after TAVR. For LCX branch in patient#1, the flowrate remained almost unchanged, however, the flow in this branch was initial significantly reduced before TAVR because of the stenosis in the middle section of the artery (peak flow for LCX was 0.062 mL/s, while for LCA and RCA is 1.5 mL/s and 0.48 mL/s respectively). For patient #2, the maximum flowrate was reduced by 19% in LAD, 17% in LCX, and 14% in RCA branches. Even after the maximum flowrate at the peak diastole, the flow reduction persists for all coronary branches during the whole diastole for both patients after TAVR (Fig. 14b). Such considerable reduction of flow could contribute to cases of ischemic lesions and promote thrombus formation.

Discussion
CAD is present in approximately 50% of the TAVR population, but this has decreased as the use of TAVR has migrated towards younger patients 83 . The question, however, of if CAD should be treated or reduced in severity prior to TAVR is still a topic for debate 83 . Additionally, TAVR and percutaneous coronary intervention (PCI) can be performed in parallel, which may reduce mortality as well as the number of vascular punctures required but may also require a larger volume of contrast agent, which could place additional strains on the kidneys 83 . Coronary arteries are supplied with blood mainly during diastole, and due to the disturbed flow associated with PVL 87 , blood entering the coronary circulation may be disrupted. The complications resulting from this is relatively unknown, and more research is needed. Hemodynamic changes, which were assessed using noninvasive computational models in this paper, may provide insight into health complications following TAVR, which may go undetected in purely anatomical examinations 88 . In the present work, there are several findings which should be individually discussed:

Improvements of coronary perfusion pressure and systemic arterial compliance after TAVR are poor indicators of coronary flow recovery in presence of paravalvular leak. AS disrupts coro-
nary flow due to the low coronary perfusion pressure 81,89 and extravascular compressive forces 81,90 , commonly associated with lower systemic arterial compliance and higher arterial resistance 81,91,92 . After TAVR, immediate increase in coronary flow is expected, as a result of increased aortic diastolic pressure (with increased forward pressure gradient at the coronary ostium) and decreased LV end diastolic pressure 81,82 . However, our findings revealed that for patients who undergo TAVR and suffer from PVL, despite the increase of aortic pressure and systemic arterial compliance, there is considerable decrease in coronary flow during diastole. We observed (Figs. 9 and 10) that in the presence of PVL, a considerable portion of the forward flow towards coronary ostium diverges towards the left ventricle, leading to a decreased coronary flow. Furthermore, in agreement with recent studies 93-95 , our results demonstrate the coronary flow is impeded if the distance between stent and coronary ostium is restricted after TAVR (Fig. 9). Such decrease in coronary blood flow is associated with reduced capacity to augment myocardial oxygenation, leading to LV dysfunction, increased apoptosis (which is linked to myocardial fibrosis and is an independent indicator of mortality) and sudden death 81,[96][97][98] .
In all patients with PVL following TAVR: No improvement of coronary flow post-TAVR. Although an increase in coronary flow is expected after AS removal and TAVR implantation 82 , our results showed that for all patients who had PVL following TAVR, a universal reduction of flow occurs during diastole for all coronary branches (Fig. 17). Recent studies suggest that despite the early improvements of systolic flow right after TAVR, coronary diastolic flow might not improve during the long-term (6-month) follow up 99 . Our results show that the coronary diastolic flow recovery is even worse for patients with PVL following TAVR. Reduced flow in coronaries could affect the outcomes of revascularization and might play a role in the pathophysiological abnormalities leading to heart failure or increased risk of cardiovascular death 96 . Sinus and neo-sinus washout after TAVR may be impaired in presence of paravalvular leak. TAVR can disturb the vortical structures inside the Valsalva sinuses, which are essential for the washout of sinus flow, assisting the smooth closure of the valve and providing flow to the coronary arteries during diastole [100][101][102] . While the sinus and neo-sinus washout efficacy of different transcatheter heart valves are still under debate 103 , our findings demonstrate that in addition to the TAVR influence on the aortic root morphology, PVL exacerbates the washout mechanism for the sinus and neo-sinus regions. We observed that the PVL jet substantially drains the flow from the sinus and neo-sinus regions, leading to pull the vortices out of the neo-sinus regions. Consequently, the vortices in the neo-sinus regions have less power to transfer the flow out of the leaflet roots. In addition, our results showed that NCC neo-sinus could be influenced the most by PVL, however, the LCC and RCC neo-sinuses irregular washout amplification depends on the severity of PVL and its location. The inefficient sinus and neo-sinus washout favors the thrombotic events after TAVR 86,100 . Subclinical leaflet thrombosis risk and hypo-attenuated leaflet thickening may be exacerbated in presence of paravalvular leak. The clinical understanding of leaflet thrombosis after TAVR www.nature.com/scientificreports/ is limited and little is known about the correlation of leaflet thrombosis with local hemodynamics 104,105 . Hypoattenuating leaflet thickening (HALT) is a thin layer of thrombus covering the aortic side of the leaflets due to subclinical leaflet thrombosis 106 . Several risk factors have been reported for thrombosis after TAVR, including reduced valve durability, restricted leaflet motion and stroke [107][108][109][110] . In addition to the agreement between our findings and previous studies 105,109,111 regarding TAVR stent morphology effect on blood stasis, we found that the PVL exacerbates the blood stasis volume in the neo-sinus regions nonuniformly and asymmetrical with respect to the valve center. While it has been reported that flow stasis risk is almost equal for LCC, RCC and NCC neo-sinuses 109 , our results revealed that not only PVL increases the blood stasis and thrombus risk in neosinus regions globally, but also is different for each neo-sinus depending on the PVL severity and location. We observed that the NCC neo-sinus is more prone to be exposed to stagnant flow and is therefore at higher risk of leaflet thrombosis than LCC and RCC.
PVL exacerbated aortic root and coronary arteries hemodynamics (local). The jets emerging from the PVL orifice substantially alters the vortical structure in the aortic root, creating disturbed flow, leading to very high shear stress at the aortic root wall. Our results demonstrate that PVL amplifies non physiological flow patterns, and consequently increases TAWSS after TAVR, especially around the leakage location. The local abnormalities in WSS are thought to stimulate aneurysm formation or lead to progressive dilation of aortic root and ascending aorta 88,112 .
On the other hand, our findings show that PVL leads to a significantly lower shear stress at the coronary walls due to the decreased blood supply during diastole after TAVR. This makes the coronary arteries susceptible to atherosclerosis, due to the low wall shear stress-induced inflammatory activation of endothelium mainly at the inner bend of curved arteries, ostia of branches and lateral walls of bifurcations 76,79 . Therefore, the decreased wall shear stress is associated with enlargement of plaque area, increased plaque eccentricity and reduced vessel area 76,[113][114][115] .
In all patients with PVL following TAVR: increased shear stress of aortic root and decreased shear stress of coronary arteries. For all patients in our study, PVL following TAVR exacerbated the shear stress during diastole (calculated through TAWSS) at aortic root and coronary arteries. TAWSS universally reduced in all branches of coronary arteries for all patients, and in contrast, significant increase of TAWSS was observed at the aortic root and around the leakage cite (Fig. 15). While the correlation of decreased shear stress at the coronaries with increased risk of plaque progression has been shown previously 76,79 , recent clinical studies also suggest that increased WSS at the aortic root could lead to ascending aorta dilation and rupture 116,117 .

PVL worsened the left ventricular hemodynamics (global).
Our results showed that moderate to severe PVL increased the burden on the LV for all patients. Despite the LV pressure reduction and increase in aortic pressure post-TAVR, LV workload increased for all patients as a result of volume overload following PVL. Therefore, PVL following the malpositioning of TAVR causes an overloaded LV, resulting in faster cardiac tissue damage and LV dysfunction. In addition, an overloaded LV may lead to other valvular diseases such as mitral regurgitation or exacerbate the existed regurgitation for patients with mixed valvular disease who receive TAVR 12,19,21 .
As shown in Fig. 16, our results showed that for all patients, the overall decrease in end diastolic pressure and increase of ascending aorta pressured lead to improved perfusion pressure. Moreover, systemic arterial compliance was improved for most of the patients (SAC reduced only for one patient (Fig. 16)). However, PVL following TAVR lead to an increased workload for most patients (LV load reduced only for one patient (Fig. 16)). The increased workload contributes to progressive myocardial fibrosis and eventually myocardial dysfunction 118,119 . Limitations of current clinical imaging modalities to capture coronary flow. Over the past decade, the use of medical imaging has drastically increased. In spite of amazing advancements in medical imaging, medical imaging on its own cannot quantify local and global hemodynamics in coronaries 120,121 . As the need for patient-specific diagnostic methods continues to be studied, understanding the strengths and limitations of imaging modalities for coronaries is critical toward creating precise diagnostic tools:(1) Computed tomography coronary angiography (CTCA): CTCA has a high spatial resolution allowing for visualization of coronary plaque and stenosis geometry 22,122 . However, CTCA suffers from temporal resolution challenges and requires the use of radiation, which is associated with health concerns especially in younger patients who need several scans throughout their lifetime 123 . CTCA does not provide any local and global hemodynamics measurements; (2) 4D flow magnetic resonance imaging (4D flow MRI): 4D flow MRI is an emerging technology to allow local hemodynamic assessment in valvular, vascular and ventricular diseases. However, use of 4D flow MRI is limited in patients with implanted medical devices as they remain a major risk during the examination. Moreover, complete and thorough analysis of local hemodynamics in coronaries is not possible 22 , due to the limited temporal resolution (4-D flow MRI has relatively high spatial resolution but lower temporal resolution (20 ms highest)). 4D flow MRI could not provide global hemodynamics; (3) Doppler echocardiography (DE): DE does not have the ability to quantify local hemodynamics through coronaries as well global hemodynamics 124 (4) Ultrafast ultrasound: Ultrafast ultrasound is an alternative option for DE, as it can image the heart at a rate of a thousand images per second 125 . Recently, ultrafast ultrasound has been combined with coronary Doppler imaging for quantification of local hemodynamics, which has aided in the diagnosis of PVR 125 . However, it has a limited imaging depth of 45 mm and cannot provide absolute quantification of flowrate for adult patients with coronary disease 125 (5)  www.nature.com/scientificreports/ none of them can provide local and global hemodynamics 128 . (6) Coronary angiography: Coronary angiography involves the transmission of a catheter into the coronary artery and the injection of a contrast medium into the bloodstream which is then viewed under X-ray examination 129 . Despite the benefits, coronary angiography is a highly invasive procedure that has shown to poorly measure FFR and evaluate the hemodynamic significance in coronaries 129 .
Limitations of current computational modeling to capture coronary flow. A clinically useful computational diagnostic framework should evaluate both global and local hemodynamics by quantifying three main requirements: (1) metrics of circulatory function (global hemodynamics), (2) metrics of cardiac function (global hemodynamics) and (3) Cardiac fluid dynamics (local hemodynamics) [18][19][20][21]33,130 . Few studies have been conducted to investigate the hemodynamic complexities after TAVR due to the presence of PVL using computational fluid dynamics (CFD) [67][68][69][131][132][133] . However, since: (1) patient-specific boundary conditions were not used; (2) hemodynamic validation was not performed; and (3) coronary arteries were excluded from the computational domain, the models developed in these studies did not satisfy the three requirements outlined in the Introduction [67][68][69][131][132][133] . In addition, several studies have recently used FSI as a promising tool for coronary arteries exclusively, since it allows consideration of the interactions of artery wall elastic behavior and blood flow mechanics, thus demonstrating its worth as a more realistic tool for numerical modelling of coronary arteries 58,[134][135][136][137][138][139][140][141] . While only a few numbers of these studies 58 coupled lumped parameter model-based boundary conditions with FSI modelling, the lumped-parameter models were not patient-specific. Moreover, all of these studies have excluded the aortic root and sinus geometry from the computational domain 58,[134][135][136][137][138][139][140] , and most of these studies have used simplified and idealized geometries of coronaries 137,138,142,143 . Exclusion of the aortic sinus at the upstream or using idealized geometry for coronaries could significantly affect the flow structure.
In this study, the requirements mentioned in the Introduction and Discussions have been examined in our study to evaluate the influence of TAVR on coronary arteries and the aortic root, when complications such as PVL or misalignment exist. In summary, our study showed that TAVR removed the aortic valve obstruction during ejection, reduced aortic valve pressure gradient and increased ejection fraction for all patients. However, considering the local flow parameters and cardiac function, all patients had adverse events after TAVR and are at high risk of heart failure. Therefore, despite the improvements of global circulatory function and clinical parameters, our results illustrating the details of local hemodynamics in these patients could partially explain how complications of TAVR could adversely increase the risk of thrombosis at aortic root and neo-sinus region of the valve leaflets, as well as plaque progression inside coronary arteries and subsequent long-term complications.

Conclusions
An optimal TAVR strategy is patient-specific, and there are varying factors that impacts the coronary hemodynamics including the global hemodynamic and circulatory system adaptation to post-TAVR environment, aortic root and aortic valve anatomical characteristics, coronary geometry, valve to coronary distance and PVL. The optimal stirring flow towards coronary arteries is diverged towards ventricle in presence of PVL and is associated with increased myocardium workload followed by progressive myocardial fibrosis and eventually myocardial dysfunction. The findings of this study suggests that exceptional consideration should be paid to the patients with paravalvular leakage after TAVR, as these patients are at higher risk of reduced coronary flow with reduced capacity to augment myocardial oxygenation, increased workload, leaflet thrombosis, plaque progression and future CAD. These complications are often asymptomatic and can lead to serious health conditions in the future and may have gone undetected if a hemodynamic assessment was not done. The scarcity of clinical trial data for complex dual pathology (CAD and AS) for the patient who undergo TAVR has urged surgeons to decide for revascularization on a case-by-case basis until further trial data. This makes the clinical endpoint and the decision for revascularization very subjective. Patient-specific computational simulations can predict the risk of post-TAVR complications such as PVL and leaflet thrombosis and their impacts on coronary hemodynamics to guide the surgeons for optimal intervention planning. The developed framework in this work is just such a tool to improve the clinical outcomes and guiding interventions for patients who receive TAVR and might be at risk of CAD over the course of time 144 .

Limitations
This study was performed and validated on 6 patients who underwent TAVR in both pre-and post-intervention states (12 cases). Future studies must consider further validation of the computational framework in a large population of AS patients in both pre-and post-intervention states, however, our results in this study demonstrate the ability of the framework to track changes in both cardiac, and vascular states. One limitation in our 3D FSI simulation is modelling only the diastole phase with TAV to be rigidly closed. It is important to note that PVL occurs only during the diastole and focusing only on diastole phase allows to simplify the simulation and reduce computational challenges and costs 67,68,71 . However, the good agreement between the FSI simulation and DE velocity data showed that this limitation does not affect the conclusions of this study. Another limitation in this study was assuming the coronary arteries fixed without the movements caused by the beating heart. However, some studies suggest that vessels dynamic motions might have negligible impacts on some parameters such as TAWSS 37 . Our computational framework is currently developed based on 6 cases and the inclusion of more cases will aid in improving the results with broader validations that could eventually be linked to patient's outcomes.

Data availability
The codes and the optimization algorithm are available from the correspondence author upon request.